#==========================================================================
# Mauerer, I. & Tutz, G. (2025). 
# "An Ordinal Item Response Model for Understanding Attitudes" 
# Sociological Methods and Research
#==========================================================================

#---------------------------------------------------------------------------------------
# This file contains the commands to run the models and create the tables and figures
# R version 4.5.1 
# Platform: x86_64-w64-mingw32/x64 (64-bit)
#----------------------------------------------------------------------------------------


#==================
# Packages
#==================
library(MultOrdRS)
library(xtable)
library(dplyr)

#==================
# Working Directory
#==================
setwd("...")

#==================
# Load Data
#==================
load("Data.RData")


#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
### Application: Attitudes toward Gender Equality (EVS)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (for details, see Supplemental Material A.1)

# SE: Sweden
# DE: Germany
# HU: Hungary


# variable townsize categorical (see Supplemental Material D)
EVS_SE$townsizecat   <-  factor(EVS_SE$townsize)
EVS_SE$townsizecat   <-  relevel(EVS_SE$townsize , ref = 1)


#--------------------
## scale covariates
#--------------------

# age
EVS_SE$age        <- scale(EVS_SE$age)
EVS_DE$age        <- scale(EVS_DE$age)
EVS_HU$age        <- scale(EVS_HU$age)
# income
EVS_SE$income     <- scale(EVS_SE$income) 
EVS_DE$income     <- scale(EVS_DE$income) 
EVS_HU$income     <- scale(EVS_HU$income) 
# townsize
EVS_SE$townsize   <- scale(EVS_SE$townsize)
EVS_DE$townsize   <- scale(EVS_DE$townsize)
EVS_HU$townsize   <- scale(EVS_HU$townsize)
# gender
EVS_SE$male       <- scale(EVS_SE$male)
EVS_DE$male       <- scale(EVS_DE$male)
EVS_HU$male       <- scale(EVS_HU$male)

#---------------------
## define formulas
#---------------------

### f.EVS_basic: without covariates
### f.EVS:       with covariates

f.EVS_basic <- as.formula(cbind(v72, v73, v74, v75, v76, v77, v78) ~ 1)
f.EVS       <- as.formula(cbind(v72, v73, v74, v75, v76, v77, v78) ~ male + age + income + townsize)

#---------------------
## fit models
#---------------------
set.seed(1860) 

## m.[COUNTRY]0: Basic Partial Credit Model
## m.[COUNTRY]1:   with Location Effects
## m.[COUNTRY]2a:  with Location & Individual Scaling
## m.[COUNTRY]2:   with Location & Scaling Effects

# SE
m.SE0  <- multordRS(f.EVS_basic, data = EVS_SE, control = ctrl.multordRS(RS = FALSE, cores = 6))
m.SE1  <- multordRS(f.EVS,       data = EVS_SE, control = ctrl.multordRS(RS = FALSE, cores = 6))
m.SE2a <- multordRS(f.EVS,       data = EVS_SE, control = ctrl.multordRS(XforRS = FALSE, cores = 6))
m.SE2  <- multordRS(f.EVS,       data = EVS_SE, control = ctrl.multordRS(cores = 6))

m.SE0
m.SE1
m.SE2a
m.SE2

# DE
m.DE0  <- multordRS(f.EVS_basic, data = EVS_DE, control = ctrl.multordRS(RS = FALSE, cores = 6))
m.DE1  <- multordRS(f.EVS,       data = EVS_DE, control = ctrl.multordRS(RS = FALSE, cores = 6))
m.DE2a <- multordRS(f.EVS,       data = EVS_DE, control = ctrl.multordRS(XforRS = FALSE, cores = 6))
m.DE2  <- multordRS(f.EVS,       data = EVS_DE, control = ctrl.multordRS(cores = 6))

m.DE0
m.DE1
m.DE2a
m.DE2

# HU
m.HU0  <- multordRS(f.EVS_basic, data = EVS_HU, control = ctrl.multordRS(RS = FALSE, cores = 6))
m.HU1  <- multordRS(f.EVS,       data = EVS_HU, control = ctrl.multordRS(RS = FALSE, cores = 6))
m.HU2a <- multordRS(f.EVS,       data = EVS_HU, control = ctrl.multordRS(XforRS = FALSE, cores = 6))
m.HU2  <- multordRS(f.EVS,       data = EVS_HU, control = ctrl.multordRS(cores = 6))

m.HU0
m.HU1
m.HU2a
m.HU2


#------------------------------------------
# Zhou (2019)'s Graded Response Model
#------------------------------------------
library(hIRT)

# SE

# survey items
y_SE <- cbind(EVS_SE$v72, EVS_SE$v73, EVS_SE$v74, EVS_SE$v75, EVS_SE$v76, EVS_SE$v77, EVS_SE$v78) 
# predictors for the mean
x_SE <- model.matrix( ~ male + age + income + townsize, EVS_SE)
# predictors for the variance
z_SE <- model.matrix( ~ male + age + income + townsize, EVS_SE)
# fitting a hierarchical graded response model
grZhou_SE <- hgrm(y_SE, x_SE, z_SE)
grZhou_SE 

# DE

# survey items
y_DE <- cbind(EVS_DE$v72, EVS_DE$v73, EVS_DE$v74, EVS_DE$v75, EVS_DE$v76, EVS_DE$v77, EVS_DE$v78) 
# predictors for the mean
x_DE <- model.matrix( ~ male + age + income + townsize, EVS_DE)
# predictors for the variance
z_DE <- model.matrix( ~ male + age + income + townsize, EVS_DE)
# fitting a hierarchical graded response model
grZhou_DE <- hgrm(y_DE, x_DE, z_DE)
grZhou_DE 

# HU

# survey items
y_HU <- cbind(EVS_HU$v72, EVS_HU$v73, EVS_HU$v74, EVS_HU$v75, EVS_HU$v76, EVS_HU$v77, EVS_HU$v78) 
# predictors for the mean
x_HU <- model.matrix( ~ male + age + income + townsize, EVS_HU)
# predictors for the variance
z_HU <- model.matrix( ~ male + age + income + townsize, EVS_HU)
# fitting a hierarchical graded response model
grZhou_HU <- hgrm(y_HU, x_HU, z_HU)
grZhou_HU 

#===============================================================================
# Table 3: Model Comparisons
#===============================================================================

#--------------------------------------------
# Table 3b: Performance of Empirical Models
#--------------------------------------------
Table3b           <- matrix(NA, 5, ncol = 4)
colnames(Table3b) <- c( "df", "Loglik.", "Loglik.", "Loglik.")
rownames(Table3b) <- c("Basic Partial Credit Model", 
                       "with Location Effects", 
                       "with Location and Individual Scaling",
                       "with Location and Scaling Effects",
                       "Zhou (2019)'s Graded Response Model")
Table3b[,1]      <- c(m.SE0$df, m.SE1$df, m.SE2a$df, m.SE2$df, 38)
Table3b[,2]      <- round(c(m.SE0$loglik, m.SE1$loglik, m.SE2a$loglik, m.SE2$loglik, grZhou_SE$log_Lik), digits=2)
Table3b[,3]      <- round(c(m.DE0$loglik, m.DE1$loglik, m.DE2a$loglik, m.DE2$loglik, grZhou_DE$log_Lik), digits=2)
Table3b[,4]      <- round(c(m.HU0$loglik, m.HU1$loglik, m.HU2a$loglik, m.HU2$loglik, grZhou_HU$log_Lik), digits=2)
Table3b

#----------------------------------
# Table 3c: Information Criterion
#----------------------------------
Table3c           <- matrix(NA, 5, ncol = 3)
colnames(Table3c) <- c("AIC", "AIC","AIC")
rownames(Table3c) <- c("Basic Partial Credit Model", 
                       "with Location Effects", 
                       "with Location and Individual Scaling",
                       "with Location and Scaling Effects",
                       "Zhou (2019)'s Graded Response Model")
# SE AIC
Table3c[,1]       <- round(c((-2 * (m.SE0$loglik)  + 2*22),  
                             (-2 * (m.SE1$loglik)  + 2*26),
                             (-2 * (m.SE2a$loglik) + 2*28),
                             (-2 * (m.SE2$loglik)  + 2*32),
                             (-2 * (grZhou_SE$log_Lik) + 2*38)), digits=2)

# DE AIC
Table3c[,2]       <- round(c((-2 * (m.DE0$loglik)  + 2*22),  
                             (-2 * (m.DE1$loglik)  + 2*26),
                             (-2 * (m.DE2a$loglik) + 2*28),
                             (-2 * (m.DE2$loglik)  + 2*32),
                             (-2 * (grZhou_DE$log_Lik) + 2*38)), digits=2)

# HU AIC
Table3c[,3]       <- round(c((-2 * (m.HU0$loglik)  + 2*22),  
                             (-2 * (m.HU1$loglik)  + 2*26),
                             (-2 * (m.HU2a$loglik) + 2*28),
                             (-2 * (m.HU2$loglik)  + 2*32),
                             (-2 * (grZhou_HU$log_Lik) + 2*38)), digits=2)
Table3c 


#----------------------------------
# Table 3d: Likelihood Ratio Tests
#-----------------------------------

# SE

## m0 vs. m1:
m.SE0_m.SE1       <- round(-2 * (m.SE0$loglik - m.SE1$loglik), digits=2)
p.val_m.SE0_m.SE1 <- round(pchisq(m.SE0_m.SE1, df = 4, lower.tail = FALSE), digits=2)

## m2a vs. m2:
m.SE2a_m.SE2       <- -2 * (m.SE2a$loglik - m.SE2$loglik)
p.val_m.SE2a_m.SE2 <- round(pchisq(m.SE2a_m.SE2, df = 4, lower.tail = FALSE), digits=2)

# DE

# m0 vs. m1:
m.DE0_m.DE1       <- -2 * (m.DE0$loglik - m.DE1$loglik)
p.val_m.DE0_m.DE1 <- round(pchisq(m.DE0_m.DE1, df = 4, lower.tail = FALSE), digits=2)

## m2a vs. m2:
m.DE2a_m.DE2       <- -2 * (m.DE2a$loglik - m.DE2$loglik)
p.val_m.DE2a_m.DE2 <- round(pchisq(m.DE2a_m.DE2, df = 4, lower.tail = FALSE), digits=2)

# HU

# m0 vs. m1:
m.HU0_m.HU1       <- -2 * (m.HU0$loglik - m.HU1$loglik)
p.val_m.HU0_m.HU1 <- round(pchisq(m.HU0_m.HU1, df = 4, lower.tail = FALSE), digits=2)

## m2a vs. m2:
m.HU2a_m.HU2       <- -2 * (m.HU2a$loglik - m.HU2$loglik)
p.val_m.HU2a_m.HU2 <- round(pchisq(m.HU2a_m.HU2, df = 4, lower.tail = FALSE), digits=2)

Table3d           <- matrix(NA, 2, ncol = 7)
colnames(Table3d) <- c("df", "chi2" ,"p-val.", "chi2" ,"p-val.", "chi2" ,"p-val.")
rownames(Table3d) <- c("gamma=0", "alpha=0")
Table3d[,1]       <- c(4,4)
Table3d[,2]       <- round(c(m.SE0_m.SE1, m.SE2a_m.SE2), digits=2)
Table3d[,3]       <- round(c(p.val_m.SE0_m.SE1, p.val_m.SE2a_m.SE2), digits=2)
Table3d[,4]       <- round(c(m.DE0_m.DE1, m.DE2a_m.DE2), digits=2)
Table3d[,5]       <- round(c(p.val_m.DE0_m.DE1, p.val_m.DE2a_m.DE2), digits=2)
Table3d[,6]       <- round(c(m.HU0_m.HU1, m.HU2a_m.HU2), digits=2)
Table3d[,7]       <- round(c(p.val_m.HU0_m.HU1, p.val_m.HU2a_m.HU2), digits=2)
Table3d

#===============================================================================
# Figure 2: Attitudes toward Gender Equality: Location and Scaling Effects
#===============================================================================

# SE
pdf(file="SE_effects.pdf", width=7, height=4.5)
par(cex.main=1.1,  mar=c(2.5,2.5,1,0), mgp=c(1.5,0.5,0))
plot(m.SE2, main = "Sweden",  alpha=0.05, ylim=c(0.7,1.75), xlim=c(0.5, 1.7), cex.lab=1, cex.axis=0.8,
     xlab="favor                                oppose",
     ylab="strong                                     weak")
dev.off()

# DE
pdf(file="DE_effects.pdf", width=7, height=4.5)
par(cex.main=1.1, mar=c(2.5,2.5,1,0), mgp=c(1.5,0.5,0))
plot(m.DE2, main = "Germany",  alpha=0.05, ylim=c(0.7,1.75), xlim=c(0.5, 1.7), cex.lab=1, cex.axis=0.8,
     xlab="favor                                 oppose",
     ylab="strong                                      weak")
dev.off()

# HU
pdf(file="HU_effects.pdf", width=7, height=4.5)
par(cex.main=1.1,  mar=c(2.5,2.5,1,0), mgp=c(1.5,0.5,0))
plot(m.HU2, main = "Hungary",  alpha=0.05, ylim=c(0.7,1.75), xlim=c(0.5, 1.7), cex.lab=1, cex.axis=0.8,
     xlab="favor                                 oppose",
     ylab="strong                                      weak")
dev.off()

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
### Application: Evaluation of Candidate Traits (ANES)
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (for details, see Supplemental Material A.2)

# democratic (dem) candidate Clinton
# republican (rep) candidate Trump

#----------------------
## scale covariates
#----------------------
ElectionUS2016$female          <- scale(ElectionUS2016$female)
ElectionUS2016$age             <- scale(ElectionUS2016$age)
ElectionUS2016$AfricanAmerican <- scale(ElectionUS2016$AfricanAmerican)
ElectionUS2016$Latino          <- scale(ElectionUS2016$Latino)
ElectionUS2016$educ            <- scale(ElectionUS2016$educ)

#---------------------
## define formulas
#---------------------

# Clinton
f.ANESdem_basic <- as.formula(cbind(leader_dem, cares_dem, know_dem, honest_dem, mind_dem, even_dem) ~ 1)
f.ANESdem       <- as.formula(cbind(leader_dem, cares_dem, know_dem, honest_dem, mind_dem, even_dem) ~ 
                                female + age + AfricanAmerican + Latino + educ)

# Trump
f.ANESrep_basic <- as.formula(cbind(leader_rep, cares_rep, know_rep, honest_rep, mind_rep, even_rep) ~ 1)
f.ANESrep       <- as.formula(cbind(leader_rep, cares_rep, know_rep, honest_rep, mind_rep, even_rep) ~ 
                                female + age + AfricanAmerican + Latino + educ)

#---------------------
## fit models
#---------------------
set.seed(1860) 
## m.[CANDIDATE]0: Basic Partial Credit Model
## m.[CANDIDATE]1:  with Location Effects
## m.[CANDIDATE]2a: with Location & Individual Scaling
## m.[CANDIDATE]2:  with Location & Scaling Effects

# Clinton
m.Clinton0  <- multordRS(f.ANESdem_basic, data = ElectionUS2016, control = ctrl.multordRS(RS = FALSE, cores = 6))
m.Clinton1  <- multordRS(f.ANESdem,       data = ElectionUS2016, control = ctrl.multordRS(RS = FALSE, cores = 6))
m.Clinton2a <- multordRS(f.ANESdem,       data = ElectionUS2016, control = ctrl.multordRS(XforRS = FALSE, cores = 6))
m.Clinton2  <- multordRS(f.ANESdem,       data = ElectionUS2016, control = ctrl.multordRS(cores = 6))

m.Clinton0
m.Clinton1
m.Clinton2a
m.Clinton2

# Trump
m.Trump0  <- multordRS(f.ANESrep_basic, data = ElectionUS2016, control = ctrl.multordRS(RS = FALSE, cores = 6))
m.Trump1  <- multordRS(f.ANESrep,       data = ElectionUS2016, control = ctrl.multordRS(RS = FALSE, cores = 6))
m.Trump2a <- multordRS(f.ANESrep,       data = ElectionUS2016, control = ctrl.multordRS(XforRS = FALSE, cores = 6))
m.Trump2  <- multordRS(f.ANESrep,       data = ElectionUS2016, control = ctrl.multordRS(cores = 6))

m.Trump0
m.Trump1
m.Trump2a
m.Trump2

#-------------------------
## Loglik.
#-------------------------

# Clinton
round(c(m.Clinton0$loglik, m.Clinton1$loglik, m.Clinton2a$loglik, m.Clinton2$loglik), digits=2)

# Trump
round(c(m.Trump0$loglik, m.Trump1$loglik, m.Trump2a$loglik, m.Trump2$loglik), digits=2)


#-------------------------
## Likelihood Ratio Tests
#-------------------------

# Clinton

## m0 vs. m1:
m.Clinton0_m.Clinton1       <- round(-2 * (m.Clinton0$loglik - m.Clinton1$loglik), digits=2)
p.val_m.Clinton0_m.Clinton1 <- round(pchisq(m.Clinton0_m.Clinton1, df = 5, lower.tail = FALSE), digits=2)

m.Clinton0_m.Clinton1
p.val_m.Clinton0_m.Clinton1

## m2a vs. m2:
m.Clinton2a_m.Clinton2       <- round(-2 * (m.Clinton2a$loglik - m.Clinton2$loglik), digits=2)
p.val_m.Clinton2a_m.Clinton2 <- round(pchisq(m.Clinton2a_m.Clinton2, df = 5, lower.tail = FALSE), digits=2)

m.Clinton2a_m.Clinton2
p.val_m.Clinton2a_m.Clinton2


# Trump

## m0 vs. m1:
m.Trump0_m.Trump1       <- round(-2 * (m.Trump0$loglik - m.Trump1$loglik), digits=2)
p.val_m.Trump0_m.Trump1 <- round(pchisq(m.Trump0_m.Trump1, df = 5, lower.tail = FALSE), digits=2)

m.Trump0_m.Trump1
p.val_m.Trump0_m.Trump1

## m2a vs. m2:
m.Trump2a_m.Trump2       <- round(-2 * (m.Trump2a$loglik - m.Trump2$loglik), digits=2)
p.val_m.Trump2a_m.Trump2 <- round(pchisq(m.Trump2a_m.Trump2, df = 5, lower.tail = FALSE), digits=2)

m.Trump2a_m.Trump2
p.val_m.Trump2a_m.Trump2


#===============================================================================
# Figure 3: Evaluations of Candidate Traits: Location and Scaling Effects
#===============================================================================

# Clinton
pdf(file="Clinton_effects.pdf", width=6, height=7)
par(cex.main=1.2,  mar=c(2.5,2.5,1,0), mgp=c(1.5,0.5,0))
plot(m.Clinton2, main = "Clinton",  alpha=0.05, ylim=c(0.7,1.17), xlim=c(0.63, 1.83),
     xlab="negative                                                     positive",
     ylab="strong                                           weak")
dev.off()


# Trump
pdf(file="Trump_effects.pdf", width=6, height=7)
par(cex.main=1.2,  mar=c(2.5,2.5,1,0), mgp=c(1.5,0.5,0))
plot(m.Trump2, main = "Trump",  alpha=0.05, ylim=c(0.7,1.17), xlim=c(0.63, 1.83), 
     xlab="negative                                                    positive",
     ylab="strong                                          weak")
dev.off()


#===============================================================================
# Table 4: Covariance and Correlation of Comprehensive Models
#===============================================================================

Table4           <- matrix(NA, 15, ncol = 4)
colnames(Table4) <- c("Covariance", "" ,"corr", "(s.e.)")
rownames(Table4) <- c("Sweden", "","", "Germany", "", "", "Hungary", "","", "Clinton", "","","Trump","","")
# SE
Table4[1,1:2]   <- round(m.SE2$Sigma[1,], digits=2)
Table4[2,1:2]   <- round(m.SE2$Sigma[2,], digits=2)
Table4[3,1]     <- round(m.SE2$se.vec[["Intercept"]], digits=2)
Table4[3,2]     <- round(m.SE2$se.vec[["RespStyle"]], digits=2)
Table4[1,3]     <- round(m.SE2$coef.vec[["Correlation"]], digits=2)
Table4[1,4]     <- round(m.SE2$se.vec[["Correlation"]], digits=2)
# DE
Table4[4,1:2]   <- round(m.DE2$Sigma[1,], digits=2)
Table4[5,1:2]   <- round(m.DE2$Sigma[2,], digits=2)
Table4[6,1]     <- round(m.DE2$se.vec[["Intercept"]], digits=2)
Table4[6,2]     <- round(m.DE2$se.vec[["RespStyle"]], digits=2)
Table4[4,3]     <- round(m.DE2$coef.vec[["Correlation"]], digits=2)
Table4[4,4]     <- round(m.DE2$se.vec[["Correlation"]], digits=2)
# HU
Table4[7,1:2]   <- round(m.HU2$Sigma[1,], digits=2)
Table4[8,1:2]   <- round(m.HU2$Sigma[2,], digits=2)
Table4[9,1]     <- round(m.HU2$se.vec[["Intercept"]], digits=2)
Table4[9,2]     <- round(m.HU2$se.vec[["RespStyle"]], digits=2)
Table4[7,3]     <- round(m.HU2$coef.vec[["Correlation"]], digits=2)
Table4[7,4]     <- round(m.HU2$se.vec[["Correlation"]], digits=2)
# Clinton
Table4[10,1:2]  <- round(m.Clinton2$Sigma[1,], digits=2)
Table4[11,1:2]  <- round(m.Clinton2$Sigma[2,], digits=2)
Table4[12,1]    <- round(m.Clinton2$se.vec[["Intercept"]], digits=2)
Table4[12,2]    <- round(m.Clinton2$se.vec[["RespStyle"]], digits=2)
Table4[10,3]    <- round(m.Clinton2$coef.vec[["Correlation"]], digits=2)
Table4[10,4]    <- round(m.Clinton2$se.vec[["Correlation"]], digits=2)
# Trump
Table4[13,1:2]  <- round(m.Trump2$Sigma[1,], digits=2)
Table4[14,1:2]  <- round(m.Trump2$Sigma[2,], digits=2)
Table4[15,1]    <- round(m.Trump2$se.vec[["Intercept"]], digits=2)
Table4[15,2]    <- round(m.Trump2$se.vec[["RespStyle"]], digits=2)
Table4[13,3]    <- round(m.Trump2$coef.vec[["Correlation"]], digits=2)
Table4[13,4]    <- round(m.Trump2$se.vec[["Correlation"]], digits=2)
Table4

#---------------------------------------------------
#### Supplementary Materials B: Estimation Tables  
#---------------------------------------------------

### B.1 Basic Partial Credit Models

#======================================================================================
# Table A1: Basic Partial Credit Model: Variance of Person Location Parameter theta_p
#=======================================================================================
TableA1           <- matrix(NA, 5, ncol = 2)
colnames(TableA1) <- c("sigmah_{theta_p}", "s.e.(sigmah_{theta_p})")
rownames(TableA1) <- c("Sweden", "Germany", "Hungary", "Clinton", "Trump")
TableA1[,1]       <- round(c(m.SE0$Sigma, m.DE0$Sigma, m.HU0$Sigma, m.Clinton0$Sigma, m.Trump0$Sigma), digits=2)
TableA1[,2]       <- round(c(m.SE0$se.vec[["Intercept"]], m.DE0$se.vec[["Intercept"]], m.HU0$se.vec[["Intercept"]],
                             m.Clinton0$se.vec[["Intercept"]], m.Trump0$se.vec[["Intercept"]] ), digits=2)
TableA1 

#=====================================================================================
# Table A2: Gender Equality Attitudes: Item Parameters of Basic Partial Credit Model 
#=====================================================================================

# Sweden
TableA2_SE           <- matrix(NA, 7, ncol = 6)
colnames(TableA2_SE) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}",
                          "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})")
rownames(TableA2_SE) <- c("v72", "v73", "v74", "v75", "v76", "v77", "v78")
TableA2_SE[,1:3]     <- round(m.SE0$beta.thresh, digits=2)
TableA2_SE[,4:6]     <- round(m.SE0$se.thresh, digits=2)
TableA2_SE 

# Germany
TableA2_DE           <- matrix(NA, 7, ncol = 6)
colnames(TableA2_DE) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}",
                          "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})")
rownames(TableA2_DE) <- c("v72", "v73", "v74", "v75", "v76", "v77", "v78")
TableA2_DE[,1:3]     <- round(m.DE0$beta.thresh, digits=2)
TableA2_DE[,4:6]     <- round(m.DE0$se.thresh, digits=2)
TableA2_DE 

# Hungary
TableA2_HU           <- matrix(NA, 7, ncol = 6)
colnames(TableA2_HU) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}",
                          "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})")
rownames(TableA2_HU) <- c("v72", "v73", "v74", "v75", "v76", "v77", "v78")
TableA2_HU[,1:3]     <- round(m.HU0$beta.thresh, digits=2)
TableA2_HU[,4:6]     <- round(m.HU0$se.thresh, digits=2)
TableA2_HU 

#===========================================================================================
# Table A3: Evaluation of Candidate Traits: Item Parameters of Basic Partial Credit Model
#===========================================================================================

# Clinton
TableA3_Clinton           <- matrix(NA, 6, ncol = 8)
colnames(TableA3_Clinton) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}", "deltah_{i4}",
                               "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})",  "s.e.(deltah_{i4})")
rownames(TableA3_Clinton) <- c("leader_dem", "cares_dem", "know_dem", "honest_dem", "mind_dem", "even_dem")
TableA3_Clinton[,1:4]     <- round(m.Clinton0$beta.thresh, digits=2)
TableA3_Clinton[,5:8]     <- round(m.Clinton0$se.thresh, digits=2)
TableA3_Clinton 

# Trump
TableA3_Trump             <- matrix(NA, 6, ncol = 8)
colnames(TableA3_Trump)   <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}", "deltah_{i4}",
                             "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})",  "s.e.(deltah_{i4})")
rownames(TableA3_Trump)   <- c("leader_rep", "cares_rep", "know_rep", "honest_rep", "mind_rep", "even_rep")
TableA3_Trump[,1:4]       <- round(m.Trump0$beta.thresh, digits=2)
TableA3_Trump[,5:8]       <- round(m.Trump0$se.thresh, digits=2)
TableA3_Trump 


### B.2 Partial Credit Models with Location and Scaling Effects


#===============================================================================
# Table A4: Effect Parameter Estimates
#===============================================================================
TableA4           <- matrix(NA, 27, ncol = 6)
colnames(TableA4) <- c("gamma","s.e." ,"gamma", "s.e.","alpha", "s.e.")
rownames(TableA4) <- c("Sweden", "male", "age", "income", "townsize", 
                       "Germany","male", "age", "income", "townsize", 
                       "Hungary","male", "age", "income", "townsize", 
                       "Clinton", "female", "age", "AfricanAmerican", "Latino", "educ",
                       "Trump", "female", "age", "AfricanAmerican", "Latino", "educ")

# SE
TableA4[2:5,1]    <- round(m.SE1$beta.X, digits=2)   
TableA4[2:5,2]    <- round(m.SE1$se.X, digits=2) 
TableA4[2:5,3]    <- round(m.SE2$beta.X, digits=2)   
TableA4[2:5,4]    <- round(m.SE2$se.X, digits=2) 
TableA4[2:5,5]    <- round(m.SE2$beta.XRS, digits=2)   
TableA4[2:5,6]    <- round(m.SE2$se.XRS, digits=2) 
# DE
TableA4[7:10,1]   <- round(m.DE1$beta.X, digits=2)   
TableA4[7:10,2]   <- round(m.DE1$se.X, digits=2) 
TableA4[7:10,3]   <- round(m.DE2$beta.X, digits=2)   
TableA4[7:10,4]   <- round(m.DE2$se.X, digits=2) 
TableA4[7:10,5]   <- round(m.DE2$beta.XRS, digits=2)   
TableA4[7:10,6]   <- round(m.DE2$se.XRS, digits=2) 
# HU
TableA4[12:15,1]  <- round(m.HU1$beta.X, digits=2)   
TableA4[12:15,2]  <- round(m.HU1$se.X, digits=2) 
TableA4[12:15,3]  <- round(m.HU2$beta.X, digits=2)   
TableA4[12:15,4]  <- round(m.HU2$se.X, digits=2) 
TableA4[12:15,5]  <- round(m.HU2$beta.XRS, digits=2)   
TableA4[12:15,6]  <- round(m.HU2$se.XRS, digits=2) 
# Clinton
TableA4[17:21,1]  <- round(m.Clinton1$beta.X, digits=2)   
TableA4[17:21,2]  <- round(m.Clinton1$se.X, digits=2) 
TableA4[17:21,3]  <- round(m.Clinton2$beta.X, digits=2)   
TableA4[17:21,4]  <- round(m.Clinton2$se.X, digits=2) 
TableA4[17:21,5]  <- round(m.Clinton2$beta.XRS, digits=2)   
TableA4[17:21,6]  <- round(m.Clinton2$se.XRS, digits=2) 
# Trump
TableA4[23:27,1]  <- round(m.Trump1$beta.X, digits=2)   
TableA4[23:27,2]  <- round(m.Trump1$se.X, digits=2) 
TableA4[23:27,3]  <- round(m.Trump2$beta.X, digits=2)   
TableA4[23:27,4]  <- round(m.Trump2$se.X, digits=2) 
TableA4[23:27,5]  <- round(m.Trump2$beta.XRS, digits=2)   
TableA4[23:27,6]  <- round(m.Trump2$se.XRS, digits=2) 
TableA4


#========================================================================================
# Table A5: Models with Location Effects: Variance of Person Location Parameter theta_p 
#========================================================================================
TableA5           <- matrix(NA, 5, ncol = 2)
colnames(TableA5) <- c("sigmah_{theta_p}^2", "s.e.")
rownames(TableA5) <- c("Sweden", "Germany", "Hungary", "Clinton", "Trump")
TableA5[,1]       <- round(c(m.SE1$Sigma, m.DE1$Sigma, m.HU1$Sigma, m.Clinton1$Sigma, m.Trump1$Sigma), digits=2)
TableA5[,2]       <- round(c(m.SE1$se.vec[["Intercept"]], m.DE1$se.vec[["Intercept"]], m.HU1$se.vec[["Intercept"]],
                       m.Clinton1$se.vec[["Intercept"]], m.Trump1$se.vec[["Intercept"]] ), digits=2)
TableA5

#======================================================================================
# Table A6: Gender Equality Attitudes: Item Parameters of Model with Location Effects
#======================================================================================

# Sweden
TableA6_SE           <- matrix(NA, 7, ncol = 6)
colnames(TableA6_SE) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}",
                          "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})")
rownames(TableA6_SE) <- c("v72", "v73", "v74", "v75", "v76", "v77", "v78")
TableA6_SE[,1:3]     <- round(m.SE1$beta.thresh, digits=2)
TableA6_SE[,4:6]     <- round(m.SE1$se.thresh, digits=2)
TableA6_SE 

# Germany
TableA6_DE           <- matrix(NA, 7, ncol = 6)
colnames(TableA6_DE) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}",
                          "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})")
rownames(TableA6_DE) <- c("v72", "v73", "v74", "v75", "v76", "v77", "v78")
TableA6_DE[,1:3]     <- round(m.DE1$beta.thresh, digits=2)
TableA6_DE[,4:6]     <- round(m.DE1$se.thresh, digits=2)
TableA6_DE

# Hungary
TableA6_HU           <- matrix(NA, 7, ncol = 6)
colnames(TableA6_HU) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}",
                          "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})")
rownames(TableA6_HU) <- c("v72", "v73", "v74", "v75", "v76", "v77", "v78")
TableA6_HU[,1:3]     <- round(m.HU1$beta.thresh, digits=2)
TableA6_HU[,4:6]     <- round(m.HU1$se.thresh, digits=2)
TableA6_HU

#=================================================================================================
# Table A7: Gender Equality Attitudes: Item Parameters of Model with Location and Scaling Effects
#=================================================================================================

# Sweden
TableA7_SE           <- matrix(NA, 7, ncol = 6)
colnames(TableA7_SE) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}",
                          "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})")
rownames(TableA7_SE) <- c("v72", "v73", "v74", "v75", "v76", "v77", "v78")
TableA7_SE[,1:3]     <- round(m.SE2$beta.thresh, digits=2)
TableA7_SE[,4:6]     <- round(m.SE2$se.thresh, digits=2)
TableA7_SE 

# Germany
TableA7_DE           <- matrix(NA, 7, ncol = 6)
colnames(TableA7_DE) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}",
                        "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})")
rownames(TableA7_DE) <- c("v72", "v73", "v74", "v75", "v76", "v77", "v78")
TableA7_DE[,1:3]     <- round(m.DE2$beta.thresh, digits=2)
TableA7_DE[,4:6]     <- round(m.DE2$se.thresh, digits=2)
TableA7_DE

# Hungary
TableA7_HU           <- matrix(NA, 7, ncol = 6)
colnames(TableA7_HU) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}",
                       "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})")
rownames(TableA7_HU) <- c("v72", "v73", "v74", "v75", "v76", "v77", "v78")
TableA7_HU[,1:3]     <- round(m.HU2$beta.thresh, digits=2)
TableA7_HU[,4:6]     <- round(m.HU2$se.thresh, digits=2)
TableA7_HU

#================================================================================
# Table A8: Candidate Traits: Item Parameters of Model with Location Effects
#================================================================================

# Clinton 
TableA8_Clinton           <- matrix(NA, 6, ncol = 8)
colnames(TableA8_Clinton) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}", "deltah_{i4}",
                               "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})",  "s.e.(deltah_{i4})")
rownames(TableA8_Clinton) <- c("leader_dem", "cares_dem", "know_dem", "honest_dem", "mind_dem", "even_dem")
TableA8_Clinton[,1:4]     <- round(m.Clinton1$beta.thresh, digits=2)
TableA8_Clinton[,5:8]     <- round(m.Clinton1$se.thresh, digits=2)
TableA8_Clinton

# Trump
TableA8_Trump             <- matrix(NA, 6, ncol = 8)
colnames(TableA8_Trump)   <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}", "deltah_{i4}",
                             "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})",  "s.e.(deltah_{i4})")
rownames(TableA8_Trump)   <- c("leader_rep", "cares_rep", "know_rep", "honest_rep", "mind_rep", "even_rep")
TableA8_Trump[,1:4]       <- round(m.Trump1$beta.thresh, digits=2)
TableA8_Trump[,5:8]       <- round(m.Trump1$se.thresh, digits=2)
TableA8_Trump


#=========================================================================================
# Table A9: Candidate Traits: Item Parameters of Model with Location and Scaling Effects
#=========================================================================================

# Clinton
TableA9_Clinton <-matrix(NA, 6, ncol = 8)
colnames(TableA9_Clinton) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}", "deltah_{i4}",
                        "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})",  "s.e.(deltah_{i4})")
rownames(TableA9_Clinton) <- c("leader_dem", "cares_dem", "know_dem", "honest_dem", "mind_dem", "even_dem")
TableA9_Clinton[,1:4] <- round(m.Clinton2$beta.thresh, digits=2)
TableA9_Clinton[,5:8] <- round(m.Clinton2$se.thresh, digits=2)
TableA9_Clinton

# Trump
TableA9_Trump           <- matrix(NA, 6, ncol = 8)
colnames(TableA9_Trump) <- c("deltah_{i1}","deltah_{i2}", "deltah_{i3}", "deltah_{i4}",
                             "s.e.(deltah_{i1})", "s.e.(deltah_{i2})", "s.e.(deltah_{i3})",  "s.e.(deltah_{i4})")
rownames(TableA9_Trump) <- c("leader_rep", "cares_rep", "know_rep", "honest_rep", "mind_rep", "even_rep")
TableA9_Trump[,1:4]     <- round(m.Trump2$beta.thresh, digits=2)
TableA9_Trump[,5:8]     <- round(m.Trump2$se.thresh, digits=2)
TableA9_Trump


### D Alternative Codings

#---------------------
## define formula
#---------------------
f.EVScat       <- as.formula(cbind(v72, v73, v74, v75, v76, v77, v78) ~ male + age + income + townsizecat)

#---------------------
## fit model
#---------------------
set.seed(1860) 

m.SE2cat <- multordRS(f.EVScat, data = EVS_SE, control = ctrl.multordRS(cores = 6))
m.SE2cat 

# Calculation of AIC
-2* (-4180.307) + 2*(38) # = 8436.614

# Model reported in the MS
-2* (-4183.11) + 2*(32) # = 8430.22


#===============================================================================
# Table A10: Effect Parameter Estimates: Townsize as Categorical (ref: 1)
#===============================================================================
TableA10           <- matrix(NA, 8, ncol = 4)
colnames(TableA10) <- c("gamma", "s.e.","alpha", "s.e.")
rownames(TableA10) <- c("Sweden", "male", "age", "income", 
                        "townsize2", "townsize3", "townsize4", "townsize5")

# SE
TableA10[2:8,1]    <- round(m.SE2cat$beta.X, digits=2)   
TableA10[2:8,2]    <- round(m.SE2cat$se.X, digits=2) 
TableA10[2:8,3]    <- round(m.SE2cat$beta.XRS, digits=2)   
TableA10[2:8,4]    <- round(m.SE2cat$se.XRS, digits=2) 
TableA10
